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ABSTRACT 


This  document  contains  information  on  the  research  accomplished  under 
AFOSR  Grant  No.  AFOSR  77-3219  during  the  time  period  1  February  1981 
through  31  January  1982.  The  work  covers  several  different  areas  of 
optical  computing  research.  Studies  of  the  limitations  of  incoherent 
optical  iterative  processors  for  inversion  of  simultaneous  linear 
equations  are  reported.  Studies  of  methods  for  finding  the  eigenvalues 
and  inverting  circulant  matrices  are  described.  Development  of  ideas 
relating  to  the  possible  use  of  optics  for  interconnections  in 
integrated  circuits  are  reported.  Finally,  new  ideas  in  the  area  of 
fiber-optic  signal  processing  are  reported.  Publications  supported  by 
the  grant  during  the  past  year  are  also  detailed. 


INTRODUCTION 


This  report  covers  the  work  performed  on  AFOSR  Grant  No.  77-3219 
during  the  time  period  1  February  1981  through  31  January  1982.  It  is 
divided  into  7  sections,  the  first  of  which  is  this  Introduction. 
Immediately  following,  we  summarize  the  results  of  a  study  of  limita¬ 
tions  of  incoherent  iterative  optical  processors,  with  the  details  pre¬ 
sented  in  an  appendix.  Section  III  contains  a  summary  on  the  progress 
on  our  continuing  studies  of  coherent  optical  methods  for  finding  the 
eigenvalues  of  circulant  matrices  and  for  inverting  those  matrices. 
Section  IV  reports  the  current  status  of  our  ideas  related  to  optical 
interconnections  for  integrated  circuits,  including  some  thoughts  on 
problems  to  which  such  techniques  might  be  applied  immediately.  Sec¬ 
tion  V  reports  the  results  of  a  short  study  of  the  multiplicative  na¬ 
ture  of  speckle  noise,  which  is  part  of  our  continuing  interest  in  the 
problem  of  suppressing  speckle  noise  in  coherently  formed  images.  Sec¬ 
tion  VI  reports  ideas  developed  over  the  past  year  on  new  approaches  to 
fiber-optic  signal  processing.  Finally,  section  VII  details  the  public¬ 
ations  supported  by  the  grant  during  the  past  year. 
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II.  LIMITATIONS  OF  ITERATIVE  INCOHERENT  OPTICAL  PROCESSORS 


A  major  effort  during  the  past  year  was  devoted  to  studies  of  the 

limitations  of  incoherent  optical  iterative  processors.  This  work  has 

been  motivated  in  large  measure  by  the  earlier  studies  of  Psaltis  and 

Casasent\  who  proposed  the  use  of  an  iterative  incoherent 

processor  for  solving  sets  of  simultaneous  linear  equations.  The 

optical  processor  is  the  incoherent  matrix- vector  multiplier  of  the 

2 

type  developed  at  Stanford  under  this  grant  ,  but  with  the  output 
fed  back  to  the  input. 

The  use  of  optical  processors  in  an  iterative  mode  raises  serious 
questions  concerning  the  effects  of  noise,  non-linearities,  and  other 
defects  that  always  plague  analog  systems.  Our  goal  was  to  analyze  the 
consequences  of  some  of  these  defects  on  the  performance  of  systems  of 
this  type.  The  results  of  this  study  are  presented  in  detail  in 
Appendix  A,  which  is  a  reprint  of  a  paper  recently  published  in 
Applied  Optics.  Thv.  problems  studied  were:  (l)  convergence 

restrictions  on  the  eigenvalues  of  the  matrix  of  coefficients  of  the 
set  of  linear  equations,  and  methods  for  relieving  these  restrictions; 

(2)  effects  of  gain  imbalance  in  the  feedforward  and  feedback  loops; 

and  (3)  effects  of  noise  introduced  by  the  successive  detection 
processes  on  the  convergence  of  the  algorithm.  Of  these  various 

problems,  by  far  the  least  explored  is  the  last. 


In  our  view,  the  most  important  of  the  results  presented  is  an 
expression  that  predicts  the  effects  of  noise  on  the  performance  of 
such  a  processor.  Suppose  that  the  equations  to  be  solved  are 
described  by  the  matrix  relation 

B  x  =  c 

where  B  is  a  matrix  of  coefficients,  x  is  a  vector  of  unknowns,  and 
c  is  a  known  vector.  Our  results  show  that  the  signal- to-noise  ratio 
at  the  output  of  the  processor  on  the  n^*1  iteration  is  bounded  by 

(S/N)n  =  (  |1  B_1c|l  )/(  )J  Kj  ) 

where  B  ^  is  the  inverse  of  the  matrix  B,  K  is  the 

n 

covariance  matrix  of  the  output  noise,  and  the  |j  signs  indicate  a 
matrix  norm.  The  norm  of  the  covariance  of  the  output  noise  was  shown 
to  be  giver  by 

i)xJI  “<s*'2  (1  -  ll?||2n+2)/0  -  )|f|(2) 

where  F  =  I  -  B,  I  being  the  identity  matrix.  The  signal  to 
noise  ratio  so  defined  can  be  shown  to  decrease  on  each  iteration, 
ultimately  approaching  a  limiting  value  after  many  iterations. 

More  work  needs  to  be  done  to  understand  the  practical  implications 
of  this  theoretical  result.  However,  we  anticipate  that  more  detailed 


examination  will  show  that  the  limiting  value  of  the  signal-to-noise 
ratio  will  depend  strongly  on  the  spread  of  the  eigenvalues  of  the 
matrix  B,  as  well  as  on  how  close  to  the  allowable  convergence 
boundaries  those  eigenvalues  lie.  More  work  is  planned  in  this  area. 

The  reader  is  referred  to  Appendix  A  for  the  full  details  of  the 
analytical  treatment. 


III.  OPTICAL  METHODS  FOR  FINDING  EIGENVALUES  OF  CIRCULANT  MATRICES 

For  the  majority  of  this  year  we  have  had  in  progress  a  research 
effort  aimed  at  developing  coherent  optical  techniques  for  finding  the 
eigenvalues  of  circulant  matrices.  The  methods  are  extendable  to  the 
problem  of  inverting  such  matrices  as  well.  We  describe  these  ideas 
here,  and  report  our  progress  over  the  past  12  months. 

The  methods  rest  on  the  well-known  fact  that  the  discrete  Fourier 
transform  (DFT)  is  the  linear  transformation  that  diagonalizes  any 
circulant  matrix.  The  complex  elements  of  the  diagonalized  matrix  are 
the  eigenvalues  of  the  original  matrix.  The  problem  then  reduces  to 
one  of  modifying  a  coherent  optical  system  such  that  the  continuous 
Fourier  transform  relations  ordinarily  obtained  between  focal  planes  of 
a  positive  lens  becomes  a  discrete  Fourier  transform.  We  have  shown  in 
another  publication^  that  such  a  modification  is  possible  if  the 
matrix  to  be  diagonalized  is  replicated  several  times  in  the  front 
focal  plane  of  a  lens,  and  if  the  light  distribution  in  the  rear  focal 
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plane  is  sampled  by  a  discrete  array  of  detectors.  The  outputs  of 
these  detectors  are  proportional  to  the  squared  magnitudes  of  the 
eigenvalues  in  question.  If  the  full  complex  values  are  desired, 
interferometric  detection  must  be  performed  in  the  rear  focal  plane. 

Once  a  distribution  of  light  amplitude  proportional  to  the 
eigenvalues  in  question  is  produced,  then  the  possibility  of  inverting 
the  matrix  arises.  If  a  coherent  optical  light  valve  with  an 
appropriate  non-linear  characteristic  (amplitude  transmittance 
proportional  to  the  logarithm  of  incident  intensity)  is  available,  then 
it  is  possible  to  show  that  the  complex  amplitude  of  the  light 
transmitted  by  such  a  light  valve  is  proportional  the  reciprocal  of  the 
original  complex  eigenvalues.  One  further  discrete  Fourier  transform 
then  produces  in  the  output  plane  a  series  of  spots,  each  with  an 
amplitude  proportional  to  one  element  of  the  inverse  matrix. 

The  conversion  of  a  coherent  optical  system  from  one  that  performs 
continuous  transforms  to  one  that  performs  discrete  transforms  requires 
replication  of  the  input  matrix,  as  mentioned  above.  In  addition,  in 

order  to  separate  the  discrete  spots  representing  the  eigenvalues  by 

* 

amounts  that  make  their  observation  easy  with  the  naked  eye,  rather 
high  resolution  and  minified  input  matrices  are  needed.  Accordingly,  a 
major  effort  was  mounted  to  develop  the  capability  of  accurately 
writing  replicated  matrices  to  photographic  film.  In  order  to 
accomplish  this  goal,  permission  was  obtained  to  use  the  DICOMET  laser 
printer  at  NASA  Ames  Research  Center.  This  printer  can  be  accessed 
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remotely  by  a  telephone- connected  computer  terminal.  The  end  result  is 
a  transparency  suitable  for  use  in  a  coherent  optical  system. 

Several  circulant  matrices  have  been  created  in  photographic  form. 
Each  element  of  the  matrix  is  encoded  as  an  amplitude  transmittance  of 
a  cell  in  an  array  of  cells.  Unless  some  form  of  holographic  encoding 
is  used,  only  non-negative  and  real  elements  can  be  allowed  in  the 
circulant  matrix.  The  matrix  is  replicated  several  times,  as  indicated 
above. 

Current  efforts  are  aimed  at  determining  the  accuracy  with  which 
the  desired  circulant  matrix  has  been  created.  To  this  end,  the 
photographic  recording  of  the  matrix  is  scanned  on  a  PDS  digital 
scanner  available  at  Stanford.  The  values  of  density  and  transmittance 
are  thus  read  from  the  transparency,  and  compared  with  the  values  that 
the  DICOMED  printer  was  instructed  to  create.  In  this  way  we  close  the 
loop  around  the  transparency  creation  process,  allowing  us  to 
understand  the  accuracy  limitations  of  the  process.  Repeatability  of 
the  exposure  values  has  been  found  to  range  from  within  4 %  to  within 
10$,  By  restricting  the  range  of  exposures  used,  we.  anticipate  being 
able  to  operate  in  the  region  of  repeatability  to  within  4^. 

Plans  for  the  future  are  centered  on  making  measurements  of  the 
eigenvalues  of  a  circulant  matrix,  and  assessing  the  error  sources  and 
error  magnitudes.  Attention  will  also  be  turned  in  the  future  to 


studies  of  the  matrix  inversion  nrocess  alluded  to  above. 


IV.  OPTICAL  INTERCONNECTIONS  IN  INTEGRATED  CIRCUITS 


Our  long-standing  interest  in  the  possible  use  of  optics  for  making 
interconnections,  either  from  chip  to  chip  or  within  a  single  chip,  has 
continued  during  the  past  year.  Vhile  our  efforts  to  interest  a 
student  member  of  the  Integrated  Circuits  Laboratory  in  working  in  this 
area  have  not  yet  been  successful,  nonetheless  some  important 
conceptual  breakthroughs  have  been  made  in  the  past  year.  Most 
important  amoung  these  has  been  the  realization  that  there  is  a  class 
of  applications  for  which  it  is  not  necessary  to  have  on  the  chip 
either  optical  sources  or  optical  modulators.  For  these  applications, 
we  need  only  have  detectors  on  the  chip,  a  task  that  seems  far  simpler 
than  the  integration  on  silicon  of  sources  or  modulators  (the  latter 
task  being  easier  with  GaAs.  The  inclusion  of  detectors  as  part  of  a 
silicon  chip  seems  straightforward  in  principle,  although  no  doubt 
practical  problems  will  be  discovered  when  such  integration  is  actually 
attempted. 

Under  what  conditions  is  it  useful  to  contemplate  placing  only 
detectors,  rather  than  both  detectors  and  sources/modulators,  on  a 
chip?  We  believe  that  such  an  approach  makes  sense  when  the  chip  must 
receive  vast  amounts  of  data  from  the  outside  world,  but  need  only 
output  small  amounts  of  data.  One  important  problem  of  this  class  has 
been  identified  during  the  past  year.  We  are  thinking  in  particular  of 
the  electronic  systolic  array,  which  is  currently  of  considerable 


interest  in  the  signal  processing  and  VLSI  communities.  Consider  a 


simple  systolic  array  designed  to  multiply  a  length  N  vector  x  by  an 
NxN  matrix  B,  yielding  a  length  N  output  vector  y.  It  is  necessary 
to  input  to  the  processor  chip  the  following  items  in  the  sequence 
described:  (l)  the  N  elements  of  the  input  vector,  sequentially  in  time 
and  on  a  single  input  channel;  (2)  the  N  elements  of  the  matrix, 
with  as  many  as  2N-1  parallel  channels,  each  carrying  in  time  sequence 
the  elements  along  one  subdiagonal  of  the  matrix.  Thus  the  data  input 
requirements  are  dominated  by  the  necessity  to  have  2N  parallel  input 
channels  for  the  matrix  elements  and  for  the  input  vector  elements;  as 
a  consequence  any  such  chip  must  have  a  large  number  of  pins  for 
inputing  data,  unless  the  chip  is  slowed  down  to  allow  some 
multiplexing  on  pins. 

For  the  same  chip,  the  data  output  requirements  are  rather  modest. 
The  N  elements  of  the  output  vector  appear  sequentially  in  time  and  can 
be  sent  to  the  outside  world  by  means  of  a  single  output  pin.  Thus  che 
number  of  pins  connecting  the  chip  to  the  outside  world  is  determined 
primarily  by  the  requirements  for  entering  data  onto  the  chip. 

Our  proposed  solution  to  this  problem  is  illustrated  in  Fig.  1. 
The  pin  connections  to  the  outside  world,  together  with  the  associated 
bonding  pads,  are  nearly  all  eliminated  by  the  use  of  optical  input 
channels.  Of  course  pins  remain  for  the  output  of  data,  for  power  to 
the  chip,  and  for  other  necessary  functions,  but  the  brunt  of  the  pin 
connection  problem  has  been  transferred  to  a  problem  of  connecting  to 
an  integrated  set  of  detectors  on  the  chip  via  a  series  of  parallel 
optical  input  channels. 
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OPTICAL  INPUT  FOR 
MATRIX  ELEMENTS 


Figure  1.  Input  of  data  to  a  systolic  array  by  means  of  optical 

connections. 

One  may  well  ask  exactly  what  has  been  gained  by  converting  from 
electronic  pin  connections  to  optical  connections.  The  problems  with 
the  optical  approach  will  not  be  fully  appreciated  until  it  has  been 
attempted  experimentally.  However,  we  can  say  the  following.  The 
problem  of  soldering  the  pin  connections  to  bonding  pads  has  been 
eliminated  for  those  pins  that  are  no  longer  necessary.  This  soldering 
process  is  said  to  be  one  of  the  most  risky  parts  of  chip  manufacture. 
We  have  not  eliminated  the  need  for  parallel  electronic  channels 
carrying  data,  for  the  optical  sources  must  be  driven  by  parallel 
electronic  channels;  rather  we  have  transferred  the  requirement  from  on 
chip  to  off  chip.  The  connections  to  the  chip  are  non-contacting,  and 
therefore  less  likely  to  fail.  The  parallel  electronics  for  inputing 
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data  can  now  be  made  with  a  larger  geometrical  scale  than  would  be 
allowed  in  the  on-chip  environment.  This  fact  brings  both  advantages 
and  disadvantages.  On  the  one  hand  the  sii,e  of  the  opto- electronic 
input  device  may  be  far  larger  than  the  chip  itself  and  more  difficult 
to  package  in  an  extremely  small  space.  On  the  other  hand,  the 
generation  of  the  connecting  optical  signals  on  a  larger  geometrical 
scale,  followed  by  a  purely  optical  demagnification,  may  make  relative 
alignment  of  the  sources  easier  than  would  be  the  case  if  all  work  were 
done  with  the  geometrical  scales  of  the  on-chip  environment. 

The  optical  approach  to  inputing  data  is  certainly  not  without  its 
own  problems.  First,  and  foremost,  there  will  be  an  alignment  problem 
of  severe  magnitude  in  attempting  to  cast  the  proper  source  images  onto 
the  proper  detectors.  Second,  the  scattering  of  light  incident  on  the 
detectors  may  be  of  sufficient  magnitude  to  lead  to  undesired  cross¬ 
talk.  Third,  it  is  uncertain  how  small  the  detectors  on  the  chip  can 
be  made.  Since  the  optical  signals  involved  can  be  of  fairly  high 
power,  there  is  reason  to  hope  that  very  small  integrated  detectors  can 
be  used. 

Currently  we  are  searching,  together  with  Prof.  James  Meindl  of  the 
Integrated  Circuits  Laboratory,  for  a  student  who  might  undertake  the 
task  of  designing  and  building  integrated  circuits  that  incorporate 
detectors  for  communicating  with  the  outside  world. 
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V.  STUDIES  OF  THE  MULTIPLICATIVE  CHARACTER  OF  SPECKLE 


The  subject  of  speckle  and  its  removal  from  coherently  formed 
images  is  one  that  has  been  of  long-standing  interest  to  us.  We  have 
spent  some  time  during  the  past  year  investigating  two  subjects  related 
to  speckle.  The  first  has  consisted  of  a  critical  study  of  the  nature 
of  speckle  noise,  and  an  examination  of  the  question  of  whether  it  is 
or  is  not  multiplicative.  A  number  of  processing  techniques  have  been 
proposed  in  the  past  that  are  based  on  the  assumption  that  a  speckled 
image  can  be  regarded  as  the  product  of  the  ideal  image  intensity 
distribution  that  would  be  obtained  if  the  object  were  not  diffuse  in 
the  fine  structure  of  its  spatial  transmission  or  reflection 
properties,  times  a  speckle  pattern  that  would  be  obtained  if  the 
macroscopic  ideal  image  were  spatially  constant  in  transmission  or 
reflection.  Similar  assumptions  have  been  frequently  made  in  modeling 
of  speckle  by  numerical  simulation.  Our  results,  which  are  attached  in 
preprint  form  as  Appendix  B,  demonstrate  that  such  a  -  model  is  valid 
only  if  the  structure  of  the  ideal  (non-diffuse)  object  is  completely 
resolvable  by  the  imaging  system  of  concern.  When  this  is  not  the  case 
(and  it  is  seldom  true,  since  most  objects  of  interest  do  have  detail 
finer  than  the  resolution  capability  of  the  optical  system  being  used), 
the  multiplicative  model  can  be  seriously  in  error.  We  refer  the 
reader  to  Appendix  B,  which  is  scheduled  for  publication  in  APPLIED 
OPTICS  in  April,  for  further  details. 


We  have  also  begun  investigation  of  some  new  digital  processing 
methods  that  we  feel  may  have  potential  for  suppressing  speckle  in 
coherently  formed  images.  The  suppression  of  speckle  is  a  highly 
desirable  goal,  for  such  effects  severely  limit  the  resolution 
obtainable  from  synthetic  aperture  radar  imagery,  from  acoustical 
imaging  systems,  and  indeed  from  any  coherent  method  of  forming  images. 
We  are  not  yet  ready  to  discuss  our  methods  in  detail  ,  but  we 
anticipate  that  by  the  time  of  next  year's  report  we  should  have  some 
interesting  results  to  describe. 


VI.  FIBER  OPTIC  SIGNAL  PROCESSING 

Our  studies  of  systolic  architectures  in  conjunction  with  the 
optical  interconnection  project  has  led  to  some  novel  concepts  in  the 
area  of  fiber-optic  signal  processing.  This  work  was  carried  out  in 
conjunction  with  Prof.  John  Shaw  of  the  Ginzton  Laboratory  at  Stanford. 
In  fact  the  experimental  evaluation  of  the  ideas  is  being  pursued  in 
Prof.  Shaw's  group,  since  he  has  the  technology  to  bring  the  ideas  to 
fruition. 

The  conceptual  advance  that  has  taken  place  in  recent  months  is  the 
origination  of  a  signal  processing  architecture  that  can  be  implemented 
using  fiber  optics  and  that  is  more  general  than  a  simple  tapped  delay 
line.  The  new  processor  is  illustrated  diagramatically  in  Fig.  2.  We 
refer  to  it  as  the  "double  helix  processor." 


The  structure  consists  of  two  single-mode  optical  fibers  that  are 
connected  together  through  a  multitude  of  optical  couplers.  Input 
signals  travel  in  one  direction  into  one  of  the  fibers,  while  output 
signals  propagate  in  the  opposite  direction  in  the  second  fiber. 
Signals  are  coupled  into  the  input  fiber  via  intensity  modulation  of  a 
laser.  The  output  intensity  modulations  are  detected  by  means  of  a 
single  detector  at  the  end  of  the  output  fiber. 

This  processor  can  be  used  in  either  a  discrete  or  a  continuous 
mode.  In  a  discrete  mode  of  operation,  the  input  laser  is  modulated  to 
emit  pulses  of  different  intensities.  A  sequence  of  N  pulses 
represents  an  input  vector  of  N  elements,  where  it  is  assumed  that  all 
elements  are  non-negative  and  real.  The  pulses  propagating  out  of  the 
second  fiber,  if  properly  windowed  in  time,  represent  the  N  components 
of  an  output  vector.  The  2N-1  coupling  coefficients  determine  the 
structure  of  the  filter  matrix  that  is  realized.  When  the  coupling 


coefficients  are  constant  in  time,  as  would  be  the  case  with  present 
fiber-optic  technology,  the  matrix  realized  is  Toeplitz,  meaning  that 
the  filter  is  time  invariant. 


It  is  also  possible 

to 

consider  this 

processor  to 

operate 

on 

continuous  functions 

of 

time,  although 

its  impulse 

response 

is 

necessarily  discrete. 

The 

continuously  modulated  input 

laser 

then 

enters  data  into  the  system,  and  the  detector  measures  a 
continuous- time  signal  at  the  output. 

There  are  a  number  of  interesting  properties  of  this  processor. 
First,  it  can  be  regarded  as  a  fiber-optic  implementation  of  a  lattice 
filter.  The  structure  incorporates  feedback,  and  as  a  consequence 
produces  an  impulse  response  that  is  theoretically  of  infinite 
duration.  Such  a  filter  can  have  both  poles  and  zeros  within  the  unit 
circle  in  the  Z-transforn  domain.  This  result  should  be  contrasted 
with  the  case  of  a  simple  tapped  delay  line,  which  produces  an  impulse 
response  of  only  finite  duration,  and  which  can  have  only  zeros  within 
the  unit  circle  in  the  Z-transform  domain.  Because  of  the  "infinite 
impulse  response"  (HR)  characteristic  of  the  double  helix  processor, 
it  should  allow  the  construction  of  much  higher-Q  filters  (for  a  given 
number  of  couplers)  than  the  "finite  impulse  response”  (FIR)  filter 
realized  by  a  simple  tapped  delay  line. 

It  is  also  worth  mentioning  that,  when  time-changing  coupler 
technology  becomes  practical  in  the  future,  the  double  helix  processor 
can  be  used  as  a  systolic  processor,  with  the  input  vector  entered  as 
pulses  on  the  input  fiber,  and  with  matrix  elements  entered  in  proper 
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time  sequence  as  time-changing  coupling  coefficients.  In  this  case  it 
is  necessary  that  the  coupling  between  the  fibers  be  weak  in  order  to 
produce  the  desired  results. 

Experimental  realization  of  a  simple  double  helix  processor  is 
being  carried  out  under  Prof.  Shaw's  funding.  Results  of  this  work 
will  be  published  in  a  letter  to  OPTICS  LETTERS,  and  credit  for 
initiating  the  theoretical  developments  will  be  given  to  AFOSR. 


VII.  PAPERS  PUBLISHED  AND  MEETING  PRESENTATIONS 

We  summarize  in  this  section  the  various  publications  and 
presentations  made  during  the  past  year  under  support  of  the  grant. 

A.  Papers  Published 

J.W.  Goodman,  A.R.  Dias,  K.M.  Johnson,  and  D.  Peri,  "Parallel  incoherent 
optical  matrix-vector  multipliers,"  Proceedings  of  the  Workshop  on 
Optical  Signal  Processing,  Texas  Tech  University,  Lubbock,  Texas, 
116-128  (1980).  (actually  published  in  1981). 

H.J.  Caulfield,  David  Dvore,  J.W.  Goodman,  and  William  T.  Rhodes, 
Eigenvector  determination  by  noncoherent  optical  methods,"  APPLIED 


OPTICS  20,  2263-2265  (1981). 


B.  Papers  Accepted  for  Publication 


J.W.  Goodman  and  Moon  Song,  "Performance  limitations  of  an  analog  method 
for  solving  simultaneous  linear  equations,"  APPLIED  OPTICS. 

M.  Tur,  K.C.  Chin,  and  J.V.  Goodman,  "When  is  Speckle  Noise 
Multiplicative?"  APPLIED  OPTICS. 

J.W.  Goodman,  "Two  Generalizations  of  the  Fourier  Transform  in 
Signal  Processing" ,  Signal  transformations  in 
S.P.I.E.  Institute  Series. 

J.W.  Goodman,  "Architectural  development  of  optical  data  processing 
systems"  (invited),  Australian  Proc.  IEE. 

The  above  paper  will  also  appear  in  the  Proceedings  of  the  Workshop  on 
Optical  Information  Processing,  held  in  Cuernevaca,  Mexico,  in 
January,  1982. 

C.  Oral  Presentations 

J.W.  Goodman,  "Architectural  development  of  optical  data  processing 
systems",  Workshop  on  Optical  Information  Processing,  Cuernevaca, 
Mexico,  January  1982. 
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Finally,  the  contributions  of  a  large  number  of  individuals  to  the 
work  reported  here  should  be  acknowledged.  Many  contributions  to  the 
"double  helix  processor"  concept  were  made  by  M.  Tur,  a  Post-Doctoral 
Fellow  partially  supported  by  the  grant.  Work  on  the  diagonalization 
and  inversion  of  circulant  matrices  was  performed  by  Q.  Cao.  The 
theoretical  developments  pertaining  to  the  limitations  of  iterative 
processors  are  due  almost  entirely  to  Moon  Song.  Contributions  to  the 
study  of  speckle  noise  were  made  by  M.  Tur  and  K.C.  Chin. 
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Performance  limitations  of  an  analog  method  for  solving 
simultaneous  linear  equations 

Joseph  W.  Goodman  and  Moon  S.  Song 

The  limitation!  inherent  in  a  recently  proposed  analog  method  for  solving  simultaneous  linear  equations  are 
examined,  and  methods  for  overcoming  some  of  these  limitations  are  discussed.  In  its  original  form  the 
method  requires  that  all  eigenvalues  of  the  matrix  of  coefficients  lie  in  a  unit  circle  centered  on  (1,0)  in  the 
complex  plane.  Proper  scaling  of  the  matrix  and  the  data  vector  extends  this  region  to  the  entire  right  half 
of  the  complex  plane  (neglecting  the  effects  of  noise).  A  modification  of  the  algorithm  is  described  that  al¬ 
lows  the  region  to  be  further  extended  to  the  entire  complex  plane.  If  the  product  of  the  gains  in  the  forward 
and  feedback  branches  is  not  unity,  the  solution  produced  by  the  algorithm  is  shown  to  be  in  error.  Finally, 
the  effects  of  noise,  which  is  inevitably  significant  in  any  analog  realization  of  the  algorithm,  are  examined. 
Noise  is  found  to  produce  a  limiting  mean  square  error  of  the  solution,  thus  preventing  perfect  convergence 
to  the  ideal  solution  vector.  A  procedure  for  determination  of  when  to  stop  the  iteration  is  proposed. 


I.  Introduction 

In  recent  publications1-2,  an  iterative  procedure  has 
been  described  that  in  principle  allows  the  solution  of 
simultaneous  linear  equations  by  means  of  repeated 
passes  through  an  incoherent  optical  system.  The 
purpose  of  this  paper  is  to  point  out  some  limitations  of 
this  particular  algorithm  and  some  methods  for  over¬ 
coming  some  of  these  limitations. 

The  method  of  concern  is  based  on  an  incoherent 
optical  matrix-vector  multiplier  described  by  Goodman 
et  al.3  Using  such  a  system,  Psaltis  et  al. 1  constructed 
an  iterative  processor  with  feedback  that  was  designed 
to  solve  simultaneous  linear  equations,  represented  in 
matrix-vector  form  by  Bx  =  c,  where  B  is  an  N  X  N 
matrix  of  known  coefficients,  c  is  a  length  N  column 
vector  of  known  data  values,  and  x  is  a  length  N  column 
vector  of  unknowns.  Figure  1  shows  a  schematic  di¬ 
agram  of  such  a  system.  An  input  vector  representing 
a  trial  solution  enters  the  system  via  a  parallel  array  of 
light-emitting  diodes  (LEDs).  The  optical  system 
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performs  inner  products  of  this  input  vector  with  the 
row  vectors  of  the  stored  matrix  mask  M.  The  resulting 
output  vector  appears  at  the  output  of  a  parallel  array 
of  photodetectors,  where  it  is  added  electronically  to  the 
known  data  vector  c.  The  result  is  passed  back  to  the 
input  of  the  system,  where  it  replaces  the  original  trial 
solution  vector.  In  certain  conditions  this  input  vector 
converges  to  the  solution  vector  desired. 

As  will  be  demonstrated,  the  iterative  procedure  in 
question  converges  only  when  the  magnitudes  of  all  the 
eigenvalues  X,  of  the  matrix  of  coefficients  lie  within  a 
unit  circle  centered  at  the  point  (1,0)  in  the  complex  X 
plane.  In  Sec.  II  we  review  these  constraints  in  prepa¬ 
ration  for  describing'certain  methods  for  relaxing  them. 
Section  III  considers  the  question  of  scaling  the  matrix 
when  its  eigenvalues  do  not  satisfy  the  necessary  con¬ 
dition,  thus  in  principle  allowing  the  eigenvalues  to  lie 
anywhere  in  the  right  half  of  the  complex  X  plane. 
Section  IV  considers  a  modification  of  the  basic  itera¬ 
tion  matrix  that  allows  extension  of  the  procedure  to 
equations  whose  matrix  has  eigenvalues  anywhere  in  the 
complex  X  plane.  Section  V  studies  the  effects  of  a  gain 
mismatch  in  the  forward  and  feedback  paths,  demon¬ 
strating  that  such  a  mismatch  leads  to  an  error  in  the 
solution.  Finally,  Sec.  VI  considers  the  effects  of  ran¬ 
dom  noise,  which  will  always  be  present  to  some  degree 
in  such  systems,  on  the  accuracy  of  the  solution.  It  is 
shown  that  with  noise  present,  the  output  of  the  system 
fails  to  converge  after  an  infinite  number  of  iterations, 
and  the  important  problem  is  how  to  choose  the  number 
of  iterations  so  that  the  solution  obtained  is  closest  to 
the  true  solution. 


502  APPLIED  OPTICS  /  Vol.  21,  No.  3  /  1  February  1982 


DATA 

vtcron 


Tig.  1.  Schematic  diagram  of  an  iterative  optical  processing 
system. 


Fig.  2.  Allowable  region  for  eigenvalues  in  the  complex  plane. 

II.  Convergence  Conditions 

Consider  the  problem  of  solving  a  set  of  linear  equa¬ 
tions  described  by  the  equation 

Bx**c,  (1) 

where  B  is  a  known  N  X  N  matrix  of  coefficients,  c  is  a 
known  column  vector  with  N  components,  and  x  is  an 
unknown  column  vector  of  N  components,  all  of  which 
we  wish  to  find. 

A  solution  for  vector  x  can  be  found  by  means  of  an 
iterative  algorithm4 

*n+l  *  (/  -  B)xn  +  c,  (2) 

which  can  be  rewritten 

(xn+,  -  B-»C)  =  (/  -  B)(x„  -  3-ic>.  (3) 

where  we  have  assumed  that  B~l  exists. 

If  we  let yn  »  x„  —  B_1c,wemay  interpret  y„  as  the 
deviation  or  error  from  the  true  solution  B-1c.  The 
error  state  equation  can  be  written 

y*+i  ■  F  y„.  <4) 

where  F  =  /  —  B,  and  the  initial  condition  is  yo  *  *o  ~ 
B~1c.  The  convergence  condition  for  this  iterative 
procedure  is  that  y„  should  vanish  as  n  increases: 

lim  y„  *  lim  Fny 0  =  0  (51 

or  equivalently 

lira  Fn  =  0.  (6) 

The  necessary  condition  for  Eq.  (6)  to  hold  is  that  all 
the  eigenvalues  of  F  be  less  than  unity  in  magnitude, 

|X,(F)|  <  1  for  1  <  i  <  N,  (7) 

or  equivalently,  since  X ,(F)  =  1  -  X,  (B), 


|1  -  X,(B)|  <  l. 

From  Eq.  (8)  we  see  that  the  eigenvalues  of  B  must  lie 
within  a  unit  circle  in  the  complex  plane  centered  at  the 
point  (1,0)  (see  Fig.  2).  In  Sec.  Ill  we  consider  the 
problem  of  scaling  the  matrix  to  extend  the  class  of 
matrices  for  which  the  algorithm  will  converge.  In  Sec. 
IV  we  consider  an  additional  modification  that  allows 
the  matrix  to  have  eigenvalues  with  both  positive  and 
negative  real  parts. 


III.  Problem  of  Scaling 

Suppose  we  wish  to  solve  a  set  of  linear  equations 
described  by 

B*x  *  c*.  (9) 

where  the  eigenvalues  of  B*  do  not  satisfy  Eq.  (8). 
Thus  |l  —  X,(B*)|  >  1.  We  introduce  a  real  scaling 
factor  <v  so  that  |l  —  X,(aB*)|  <  1.  Now  let  B  =  aB*, 
c  =  ac*,  yielding 

(aB*)x  =  (ac*).  (10) 

This  new  equation  will  have  the  same  solution  vector 
as  the  original  Eq.  ( 1 ).  Thus  by  applying  the  iterative 
procedure  with  iteration  matrix, 

f  =  /-aB»,  (11) 

and  simultaneously  scaling  the  data  vector  by  the  same 
amount,  the  desired  solution  should  be  obtained.  By 
use  of  the  scaling  factor  a,  the  allowable  region  for  the 
eigenvalues  of  B  can  be  extended  to  the  entire  right  half 
of  the  complex  X  plane  (neglecting  restrictions  posed  by 
noise  and  speed  of  convergence). 

While  the  above  procedure  appears  straightforward, 
we  have  not  yet  addressed  the  important  question  as  to 
how  the  scaling  factor  a  should  be  chosen.  The  rate  of 
convergence  of  the  algorithm  is  dependent  on  the 
magnitude  of  the  largest  eigenvalue,  which  is  known  as 
the  spectral  radius,  of  the  iteration  matrix  F.  There¬ 
fore,  to  assure  most  rapid  convergence,  we  should  choose 
the  scaling  factor  to  make  the  magnitude  of  the  maxi¬ 
mum  eigenvalue  of  the  iteration  matrix  as  small  as 
possible.  Note  that  small  eigenvalues  of  F  are  not 
achieved  simply  by  choosing  a  small  scaling  factor,  for 
as  «  approaches  zero  all  the  eigenvalues  of  F  approach 
unity.  The  convergence  of  the  process  is  also  limited 
by  the  presence  of  random  measurement  noise.  The 
smaller  the  eigenvalues  of  F,  the  more  significant  the 
effects  of  this  noise  become.  We  defer  a  discussion  of 
noise  effects  to  Sec.  VI  and  consider  here  only  the  op¬ 
timum  choice  of  o  to  achieve  a  minimum  of  the  magni¬ 
tude  of  the  largest  eigenvalue  of  F.  At  this  point  we 
restrict  attention  to  the  class  of  Hermitian  matrices,  so 
that  all  the  eigenvalues  are  real.  For  non-Hermitian 
matrices,  we  can  use  the  transformation  method  de¬ 
scribed  in  Sec.  IV. 

Noting  that 

A,(fl)  =  A,(afl*)  “  aA,(fl*).  (12) 

we  see  that  in  passing  from  B*  to  the  iteration  matrix 
F,  the  maximum  eigenvalue  Xmax(B*)  will  be  mapped 
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into  1  “  oXro«(B*),  while  the  minimum  eigenvalue 
Xa,intB*)  will  be  mapped into  1  ~  “Xmin(B*).  Thus  to 
make  the  maximum  eigenvalue  of  Fas  small  as  possible, 
we  should  choose  a  so  that  the  maximum  and  minimum 
eigenvalues  of  F  are  equal  in  magnitude  and  opposite 
in  sign, 

(13) 

which  yields  an  optimum  value  of  a  given  by 
_ 2 _ 

“0Pt  “  Xauu(B*)  +  XmiaCB*)  '  (W> 

Now  suppose  that  the  spread  of  eigenvalues  of  B*  is 
fixed  at  value  T  =  Xmin(B*)Amax(B*)-  With  the  use  of 
the  optimum  value  of  a  above,  we  And  that  the  magni¬ 
tude  of  the  maximum  eigenvalue  of  the  iteration  matrix 
F  is  given  by 

|X«*(F)|  -  |1  -  aoptX^B*)!  -  |^|  •  (15) 

This  result  represents  the  smallest  achievable  value  of 
the  magnitude  of  the  eigenvalue  of  F  having  largest 
magnitude  and  can  be  realized  only  by  choosing  the 
scaling  constant  a  to  be  given  by  the  value  specified  in 
Eq.  (14). 

While  the  above  choice  of  a  gives  the  most  rapid 
convergence  rate  theoretically,  we  do  not  have  the  ei¬ 
genvalues  at  hand  in  advance.  But  we  can  make  an 
estimate  of  the  maximum  eigenvalue  of  B  by  using  the 
well-known  Gerschgorin  theorem5  and  the  Schur  the¬ 
orem.6  The  Gerschgorin  theorem  states 

N 

Xau  —  X bound,  i  ~  max;  1]  |  bjj  1 ,  (16a) 

i- 1 . 

N 

*max  —  Xbound,2  *  ®a *j  £  (16b) 

*-l 

The  Schur  bound  gives 

Xmu  ^  Xbound.3  “  (  2Z  2Z  5?j|  •  (16c) 

Vi-l  /-l  I 

Since  these  bounds  do  not  impose  a  heavy  computa¬ 
tional  load,  it  is  reasonable  to  compute  all  of  them  and 
choose  the  smallest  one: 

Xqui  -  X bound  “  mjn | Abound,  1 1 X bound ,2< X bound ^i-  (16d) 

But  since  it  is  not  easy  to  estimate  the  upper  bound  of 
the  minimum  eigenvalue,  it  seems  practically  reasonable 
to  choose  the  scaling  factor  as7 

a  -  ZAXbmaK,).  (17) 

IV.  Extension  to  Matrices  B  with  Eigenvalues  not 
Conflnod  to  the  Right  Half  of  the  Complex  Plane 

We  assumed  in  the  previous  discussion  that  all  the 
eigenvalues  of  the  matrix  B  lie  in  the  right  half  of  the 
complex  plane.  Here  we  discuss  further  generalizations 
that  allow  matrices  not  satisfying  this  restriction  to  be 
dealt  with. 

First  we  note  that  if  the  matrix  B  has  eigenvalues  that 
lie  within  a  unit  circle  centered  at  (—1,0)  in  the  complex 
plane,  a  simple  modification  of  the  algorithm  to 

*n+l  *  (/  +  fl)xn  +  C  (18) 


Fig.  3,  Block  diagram  of  iterative  system. 


allows  such  matrices  to  be  dealt  with.  By  introducing 
a  scaling  factor,  the  region  of  allowable  eigenvalues  can 
be  extended  to  the  entire  left  half-plane  (neglecting 
limitations  posed  by  noise  and  convergence  speed). 

Suppose  that  the  matrix  B  has  some  eigenvalues  with 
positive  real  parts  and  some  eigenvalues  with  negative 
real  parts.  Neither  algorithm  presented  will  allow  such 
matrices  to  be  included.  However,  there  is  a  simple 
modification  of  the  algorithm  that  will  do  so.  Let  the 
equation  to  be  solved,  i.e.,  Eq.  (9),  be  modified  by 
multiplying  the  left  and  right  sides  by  the  Hermitian 
transpose  B*H  of  B*,  giving 

(B*)"B*x  *  (B*)wc*.  (19) 

The  matrix  ( B*)HB *  is  non-negative  definite  (practi¬ 
cally  positive  definite  since  B*  is  assumed  nonsingular), 
and  therefore  the  use  of  an  iteration  matrix  F  =  /  — 
a(B*)HB*  and  a  proper  choice  of  a  will  assure  that  the 
eigenvalues  of  F  are  non-negative  and  less  than  unity 
in  magnitude.  Furthermore,  the  solution  set  for  Eq. 
(19)  is  precisely  the  same  as  the  solution  set  for  Eq.  (9). 
Thus  the  multiplication  of  both  the  matrix  B*  and  the 
data  vector  c*  by  (B*)H  has  generated  a  new  equation 
with  the  same  solution  set  but  with  matrix  eigenvalues 
satisfying  the  necessary  requirements.  In  this  way  the 
space  of  allowable  eigenvalues  has  been  extended  to  the 
entire  complex  plane. 

V.  Effects  of  Gain  Mismatch  in  the  Forward  and 
Feedback  Paths 

Since  the  optical  system  in  question  involves  opti- 
cal-to-electronic  and  electronic-to-optica!  conversion, 
and  since  they  are  analog  systems,  the  strong  possibility 
exists  for  differences  of  the  gains  in  the  forward  and 
feedback  paths.  Figure  3  illustrates  the  general  nature 
of  the  systems  of  concern.  G\  represents  the  inherent 
gain  matrix  associated  with  passage  of  the  signals  in  the 
forward  direction,  while  G 2  is  the  corresponding  gain 
matrix  in  the  feedback  direction.  The  matrices  are 
diagonal,  but  the  elements  along  the  diagonal  need  not 
be  equal,  since  the  different  channels  of  an  incoherent 
matrix-vector  multiplier  may  have  different  gains. 

The  presence  of  the  gain  matrices  G 1  and  G’2  modifies 
Eq.  (2)  describing  the  iterative  algorithm  to  become 

xn+i  =  GjC-,(/  —  B)  x„  +  c.  (20) 

To  see  the  effects  of  a  gain  mismatch,  suppose  that  the 
product  of  gain  matrices  is  close  to  the  identity  ma¬ 
trix, 

G2G,  «/+«.  (r”) 
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where  the  diagonal  matrix  e  represents  the  gain  mis¬ 
matches  in  the  various  channels  and  is  assumed  to  have 
small  elements.  The  iteration  equation  becomes 

x^+i  -  (/  + «)(/  -  B)x„  +  c  -  |i  -[£-<</  -  S)]|x.  +  c. 

(22) 

From  the  above  equation  we  can  see  that,  when  gain 
mismatches  exist,  the  iterative  procedure  will  produce 
it  solution  that  is  appropriate  for  a  matrix  B  -  e(I  -  B) 
rather  than  the  matrix  B.  We  now  investigate  the 
magnitude  of  the  errors  associated  with  the  solution 
vector  x. 

Suppose  that  the  iteration  procedure  has  reached 
convergence;  i.e., 

*n+l  ■**«“*.  (23) 

Then  the  iteration  equation  can  be  written 

(/  —  (/+ «)(/  —  B)]x  «  c.  (24) 

Solving  this  equation  for  x,  we  find 

x  -  (/  -  B-UU  -  BJJ-'B-'c.  (25) 

If  the  elements  of  matrix  e  are  all  small,  the  result  for  x 
can  be  approximated  as 

x  *  B-‘c  +  £->«(/  -  B)B-»c.  (26) 

Thus  in  the  presence  of  gain  mismatch,  the  iterative 
procedure  converges  to  an  incorrect  answer.  The  error 
vector  is  given  by 

Xtnw  •  B~lt(I  —  B)B_lc.  (27) 

The  relative  size  of  the  error  due  to  gain  mismatch  can 
be  written  (assuming  a  symmetric  B) 

ls„rJ  ilB-Mf-B)B->cll 

|B-‘cl  "  HB-»ci 

l|B-M/-B)11-||B-»c| 

iB->cll 

<  l|<l-llfl-‘-/||  <  u<8  -  (HB-M1  +  1),  (28) 

where  IMI  indicates  a  matrix  norm,  defined  by 

lift II  s  | largest  eigenvalue  of  KHK\ 1/2,  (29) 

and  the  superscript  H  indicates  a  Hermitian  transpose. 
Thus  as  might  be  expected,  the  bound  on  the  relative 
error  is  proportional  to  the  gain  mismatch  ||c||  and  to 
llfl-l|l  (which  is  essentially  the  inverse  of  the  smallest 
eigenvalue  of  B). 

VI.  Effects  of  Noise  on  the  Convergence  of  the 
Algorithm 

When  an  iterative  procedure  such  as  that  of  interest 
here  is  implemented  with  an  analog  system,  concern 
naturally  arises  as  to  the  effects  of  noise.  Noise  is  in¬ 
troduced  on  every  pass  through  the  loop,  and  it  should 
be  expected  that  its  effects  will  eventually  build  up  to 
the  point  where  they  cannot  be  neglected.  Our  purpose 
here  is  to  examine  the  effects  of  noise  in  as  quantitative 
a  way  as  possible. 

We  make  the  approximation  that  the  noise  involved 
is  purely  additive,  and  we  further  assume  that  it  has  zero 


W 

a 


Fig.  4.  Block  diagram  of  system  with  noise. 

mean.  Figure  4  shows  a  block  diagram  of  the  system 
considered.  The  state  equation  with  a  noise  component 
present  must  be  rewritten  as 


*n+ 1  -  (/  -  B)x„  +  c  +  w„,  (30) 

where  w„  is  the  noise  vector  of  concern.  The  error  state 
equation  becomes 

y,+i  =  Fyn  +  w„,  (31) 

where  again  we  have  used  the  definition  y„  =  x„  — 
B~le. 

The  solution  of  Eq.  (30)  can  be  divided  into  two 
components:  (1)  a  deterministic  homogeneous  solution 
S'n  due  to  the  initial  condition,  and  (2)  a  stochastic 
particular  solution  y„  due  to  the  noise  input.  Since  the 
system  is  linear,  we  may  consider  these  components 
separately.  Note  that  yn  is  the  solution  to 

9n*i  -  F$n  (32) 

with  $o  *  Xo  —  B-1c,  while  y„  is  the  solution  to 

y„+i  =  Fy„  +  w„  (33) 

with  yo  =  0.  The  complete  solution  is  given  by 

y  n=9n  +  9n  (34) 

The  homogeneous  part  of  the  solution  was  analyzed  in 
the  earlier  sections.  Here  we  consider  only  the  noisy 
part  y„. 

From  Eq.  (33)  we  write 

y»+i  =  F9n  +  (35) 

*-  o 

Squaring  both  sides  of  this  equation  and  taking  averages 
with  the  assumption  that  the  noise  is  both  white  and 
stationary  yield  the  following  expression,  for  the  co- 
variance  matrix  of  the  noise  output. 

Kn+,  ^  £(yn+1y«+1)  =  *2  i  ( FFH)\  (36) 

*-  o 

where  cr2l  is  the  covariance  of  noise  vector  w„. 

Using  the  definition  of  matrix  norm  in  Eq.  (29),  the 
following  inequality  can  be  used  to  put  a  bound  on  the 
error  covariance: 


n  1  -  |! f|i2"  +  2 

IK, I  <  *2  L  Ilf  II2*  =  *2- — 77^7"  •  (37) 

*-o  1  -  llflr 

With  the  above  result  in  mind,  we  define  a  time 
varying  (SNR)„  by 


(SNR)n  = 


IIB-'cll 

V  ilKn  il 


(33) 
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F\g.  5.  Block  diagram  of  system  with  criterion  to  stop. 


From  Eq.  (38)  we  see  that,  due  to  the  increasing  nature 
of  NfCj  with  iteration  number,  the  (SNR)„  becomes 
smaller  at  each  iteration.  Furthermore,  when  the  norm 
||F||  is  less  than  one,  the  (SNR)„  ultimately  approaches 
a  limiting  value  that  is  independent  of  the  iteration 
number. 

During  the  initial  iterations,  the  starting  condition 
will  be  sufficiently  far  from  the  true  solution  that  the 
deterministic  component  of  error  will  dominate,  and  the 
output  vector  will  start  to  move  toward  the  solution 
vector.  The  erro;  will  therefore  begin  to  decrease.  As 
the  iteration  number  grows,  the  noise  component  of  the 
output  builds  up,  and  eventually  it  becomes  the  domi¬ 
nant  component  of  error,  preventing  the  error  from 
vanishing.  Thus  the  process  initially  converges,  but 
then  the  mean  square  error  approaches  a  nonzero  as¬ 
ymptotic  value.  A  key  question  concerns  the  means  for 
judging  when  to  stop  the  iteration.  When  the  process 
is  sufficiently  slow  to  allow  a  human  observer  to  interact 
with  the  system,  judgment  can  perhaps  be  made  as  to 
when  to  terminate  the  procedure.  However,  interest 
is  greater  when  the  processing  speed  is  high,  and  in  such 
cases  the  possibility  of  human  intervention  vanishes. 
An  automatic  means  for  judging  when  to  terminate  the 
iteration  is  needed.  We  now  describe  one  such 
method. 

Suppose  we  obtain  on  a  particular  iteration  a  solution 
x  which  may  differ  from  the  true  solution  B~lc.  Ob¬ 
viously,  it.  is  our  hope  that  x  will  satisfy  the  given 
equation  as  closely  as  possible,  i.e.,  that 

BX  =  c.  (39) 

Thus  oui-  efforts  to  solve  the  equation  are  equivalent  to 
the  minimization  of  the  equation  error,  c  —  Bx,  which 
will  be  zero  if  x  happens  to  be  B_1  c.  Hence  we  can  use 
this  error  as  a  measure  of  closeness  to  the  solution.  Let 
the  equation  error  associated  with  solution  xn  be  rep¬ 
resented  by  en ;  i.e., 

e„=c-Bx„.  (40) 

Then  it  is  easy  to  demonstrate  that  e„  satisfies 

*  Fe„  (41) 

with  the  initial  condition  eo  =  c  —  Bx o-  If  B  and  c  have 
been  properly  scaled  so  that  F  is  a  stable  matrix,  we  see 
that  en  should  vanish  as  the  iteration  proceeds 

lira  e„  =  lim  Fne0  =  0.  (42) 


We  can  monitor  the  equation  error  at  each  iteration  by 
measuring  some  form  of  the  vector  norm  of  e„  (One 
easy  method  would  be  by  taking  the  absolute  value  of 
the  maximum  component  of  e„.)  This  quantity  thus 
provides  the  required  estimate  of  closeness  to  the  true 
solution. 

The  last  question  is:  how  can  we  obtain  the  equation 
error  experimentally?  If  we  note  that 

Xn+i  ■  Xn  -Bx„  +  c  »  x„  +  e„,  (43) 

which  is  readily  available  after  the  (n  +  l)st  iteration, 
we  can  obtain  en  simply  by  subtracting  x„  from  xn+1.’ 
The  system  then  must  contain  a  means  for  storing  one 
vector  x„.  A  block  diagram  of  the  method  is  illustrated 
in  Fig.  5. 

VII.  Conclusions 

We  have  discussed  some  performance  limitations  of 
a  certain  optical  implementation  of  an  iterative  algo¬ 
rithm  that  can  be  used  to  solve  simultaneous  linear 
equations.  The  method,  as  originally  proposed,  re¬ 
quires  that  all  eigenvalues  of  the  matrix  of  coefficients 
lie  within  a  circle  in  the  complex  plane  having  a  unit 
radius  and  centered  on  the  point  (1,0).  We  have  de¬ 
scribed  a  method  for  modifing  the  matrix  and  the  data 
vector  so  that  the  eigenvalues  of  the  matrix  of  interest 
can  (in  principle)  lie  anywhere  in  the  complex  plane.  It 
has  been  shown  that  a  difference  in  gains  associated 
with  the  forward  and  feedback  paths  causes  the  algo¬ 
rithm  to  converge  to  an  incorrect  solution  vector.  The 
effects  of  measurement  noise  on  the  algorithm  have  also 
been  examined.  It  has  been  shown  that  the  presence 
of  noise  affects  the  convergence  of  the  algorithm,  leading 
to  a  mean  square  error  between  the  true  solution  vector 
and  the  experimental  solution  vector  that  fails  to  vanish 
as  the  iteration  number  grows  large.  Finally,  we  have 
suggested  a  means  for  determining  when  the  iterations 
should  be  stopped,  based  on  the  difference  between  the 
solution  vectors  on  two  successive  iterations.  All  the 
above  considerations  are  of  considerable  importance 
when  the  algorithm  is  implemented  by  means  of  an 
analog  optical  system. 

The  work  described  here  was  supported  by  the  Air 
Force  Office  of  Scientific  Research. 
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APPENDIX  B 


WHEN  IS  SPECKLE  NOISE  MULTIPLICATIVE? 

M.  Tar 
K  C.  Chin * 

J.  W.  Goodman 

Information  System  Laboratory,  Stanford  University,  Stanford,  Ca.  94305 

In  coherent  illumination,  objects  with  roughness  on  the  order  of  a 
wavelength  cause  speckle  to  appear  in  their  images  as  formed  by  imaging  sys¬ 
tems  which  cannot  resolve  the  microscale  of  the  objects’  roughness  [l].  Thus, 
when  a  laser  illuminates  a  composite  object  composed  of  a  diffuser  (whose  com¬ 
plex  amplitude  transmittance  is  d( x.y'))  in  contact  with  a  transparency 
t(x',y)  (see  Fig.  1)  the  spatial  intensity  distribution  in  the  image  plane  will  be 
very  noisy  (i.e.  speckled  ).  provided  the  amplitude  point  spread  function  of  the 
optical  system  is  broad  by  comparison  with  the  microscopic  surface  variations 
of  the  diffuser.  Recent  publications  [2]-[4]  as  well  as  older  ones  [5]-[7]  have 
assumed  that  speckle  noise  is  multiplicative,  i.  e. 

^u(x.y)  =  a  4* (x.y)  /s (x.y)  (1) 

where:  Its(x ,y)  is  the  (random)  spatial  intensity  distribution  in  the  speckled 
image  of  the  transparency  £(x\y)  .  /s(x,y)  is  the  intensity  distribution  in  the 
image  of  the  diffuser  alone  (  t (x  ,y  )  =  1  )  and  /me (x.y)  is  the  incoherent  image 
of  the  transparency  t(x',y ’)  .  The  proportionality  factor  a  depends  on  the  sys¬ 
tem  parameters.  Eq.  (1)  is  mainly  based  on  the  work  of  Lowenthal  and  Arsenault 
[8],  who  showed  that  when  d(x'  ,y)  is  a  stationary,  6  correlated  random  pro¬ 
cess  with  independent  real  and  imaginary  parts  having  zero  means  and  the  same 
On  leave  from  the  Dept,  of  Physics,  Pe- f  -ig  University,  People's  Republic  of  China. 


variances,  the  mean  and  standard  deviation  of  Iu  (x  ,y  )  are  equal  to  (x  ,y  ) 
That  is. 


<Ju{x.y)>  =  4ic (z,y)  ;  <lg(x.y)>  =  2 I&c(x,y)  (2) 

where  <  >  denote  an  ensemble  average  over  different  realizations  of  the 
diffuser  and/or  the  illumination.  Using  Eqs.  (l)-(2),  we  also  find  that  a  in  Eq.  (1) 
is  given  by:  a=  <Is(x ,y)>~1  .  While  Lim  and  Nawab  [2]  specifically  restricted 
their  multiplicative  model  to  the  case  when  the  degraded  image  has  been  sam¬ 
pled  coarsely  enough  such  that  the  degradation  at  any  point  cam  be  assumed  to 
be  independent  of  that  at  all  other  points,  most  other  investigators  have  used 
Eq.  (l)  without  further  limitations. 

It  is  the  purpose  of  this  letter  to  point  out  that  Eq.  (l)  is  only  an  approxima¬ 
tion  which  is  certainly  not  valid  when  the  transparency  t{x',y)  has  spatial 
details  which  cannot  be  resolved  by  the  coherent  system.  This  observation  is  of 
practical  importance  since  most  objects  contain  fine  details  well  beyond  the 
resolution  capabilities  of  the  systems  that  are  used  to  image  them. 

Due  to  the  finite  resolving  power  of  the  imaging  system,  the  complex  wave 
amplitude  U  at  the  image  point  P  (see  Fig.  l)  is  the  result  of  a  coherent  addi¬ 
tion  of  contributions  from  many  independent  areas  within  a  finite  patch  RctU 
which  is  of  the  order  of  a  resolution  cell.  The  (random)  intensity  at  P=(xp,yp)  is 
given  by 

JtM^Xp.yp)  =  \M{Xp,yp)\z  =  |  f f  dx'  dy'  h(xp-x',yp-y')  t(x'.y’)  exp[ip(x',y')]  1(3) 

|  | 

where  h{x,y)  is  the  amplitude  impulse  response  of  the  system  (assuming  unity 

magnification)  and  <p{x’,y  )  is  a  random  phase  delay  which  is  introduced  by  the 

diffuser  at  the  object  point  (z  ,y)  ;  d(x\y  )  =  exp[i<p(x'y .  Obviously,  when 

,y)  does  not  change  appreciably  within  Rc„u  ,  Eq.  (3)  reduces  to 
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,2 

=  1*^.2/*)  1 2  ffd*  d2/'/i(xp -x'1%) V)  exp[i^(z'.v')]  (4) 

"Mil  | 

Since  in  this  Limiting  case  Imc(xp.yp)  a  \t{xp,yp)\z  and  the  second  factor  on 
the  right  hand  side  of  Eq.  (4)  is  /,(a^,T/p)  ,  Eq.  (4)  has  the  same  form  as  Eq.  (l). 
This  explains  why  Lee  [3]  was  successful  in  verifying  Eq.  (1)  for  the  flat  areas  in 
his  synthetic  aperture  radar  images.  However,  when  t{z,y)  has  spatial  details 
smaller  than  or  comparable  to  i?c#u  .  Eq.  (1)  is  no  longer  valid.  As  an  example, 
consider  the  case  of  a  sharp  opaque  edge  object  and  denote  by  /to  the  speckled 
(random)  intensity  at  the  geometrical  image  of  the  edge  G  (see  Fig.  2)  and  by 
I,  the  intensity  of  the  speckle  at  the  same  point  in  the  absence  of  the  sharp 
edge.  Based  on  a  random  walk  model  [9],  /s  is  the  intensity  of  the  sum  of  two 
independent  random  walks  Z*  and  Z~  originated  respectively  from  R£u  and 
Rchu  in  Fig.  2.  Obviously.  Iu  is  the  intensity  of  Z*  .  Assuming  exponential  dis¬ 
tributions  for  the  intensities,  it  is  readily  shown  that  the  joint  probability  den¬ 
sity  function  for /to  and  Is  is  given  by: 


pV.Jts) 


1 _ 

[  </,>-</ls>]</i9>  exp 


_  T  p _ >/■/»  Its _ 

[  </,>-</*>  ]</*>]  °[  [  </s>-</,s>  ] 

(5) 


where  <Ia  >  and  </to>  are  the  ensemble  averages  of  /,'  and  Its  ,  and  /0  is  a 
modified  Bessel  function  of  the  first  kind,  zero  order.  Therefore,  the  quotient 


Q  = 


.assumed  constant  in  Eq.  (l),  is  instead  a  continuously  distri¬ 


buted  quantity  with  a  probability  density  given  by: 


p(g)  =  zq  + 


where  B  - 


—  r~  (In  the  simplest  case,  when  RceU  is  equally  divided  by  the 


edge,  B- 2).  Since  for  large  values  of  Q,p{Q)ssQ  3,  Q  has  an  infinite 
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variance  which  means  that  the  zeroes  of  Ita  occur  independently  of  those  of 
/,  .  in  contradiction  to  Eq.  (l). 

To  further  validate  our  assertion,  we  have  carried  out  one-dimensional 
numerical  simulations  of  the  image  forming  system  of  Fig.  1  for  slanted  objects 
with  different  slopes.  The  speckled  image  was  generated  using  convolution  tech¬ 
niques  rather  than  Fourier  domain  methods.  Results  for  -p—  near  the  point  G 

in  Fig.  2,  show  that  for  two  different  realizations  of  the  diffuser,  the  steeper  the 
edge,  the  larger  the  difference  between  the  two  realizations.  We  have  also 


confirmed  that  for  B 


- 


>  assumes  the  value  1.91  which  is  the  first 


moment  of  the  distribution  p(Q)  ,  Eq.  (6). 

A  statistical  estimate  for  the  error  introduced  by  the  multiplicative  model 
can  be  obtained  from  the  mean  square  error  M[x,y)  : 

M{x,y)  =  <[/to(x,y)  -  a  lincix.y)  /,(x.y)]2>  (7) 

In  order  to  evaluate  M  we  follow  the  assumptions  and  results  of  Lowenthal  and 
Arsenault  [8],  so  that  Eq.  (7)  reduces  to: 

M/.  orZ  Jo  1 

U{X-V)  ~  -  Or.^'<£(x;y)>j  <s> 

Using  the  above  assumptions  together  with  the  additional  fact  that  the  complex 
wave  amplitudes  and  U,  at  the  image  point  (x.y)  ,  originate  from  the 
same  diffuser,  we  may  conclude  that  Uts(x,y)  and  U,(x,y)  are  two  complex 
mutually  circular  Gaussian  random  variables,  and  therefore: 


<h,(*.y) Ia{x.y)>  =  <Uis(x,y)Ute(x.y)Us(x,y)U/(x.y)> 

-  -A/ic  ix-y)  <Ia{x.y)>  +  |<Uf5 (x ,y )U/(x ,y)>  j 


I 
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The  star  denotes  complex  conjugation.  Substituting  Eq.  (9)  in  Eq.  (0).  we  find: 


M{x,y )  =  2  I&cix.y) 


1  - 


<Ut,(x.i/)Ua*(x.y)>|2 


Ane  (x,y)  <J,{x.y)> 


(10) 


Recalling  that  <d.{x\,y\)d{x'z,y2)>  =  6{x\-x‘z)6{y\^y‘z)  and  that 
Idfx'.y’)!  =  1  .  Eq.  (10)  may  be  rewritten  as: 


M(x,y)___  i  _ 
2/^c(*.y) 


1/ 


|Z 


dx  dyt{x,y')  |/i(x-x',y-y')  |2 


+mm  +«• 

f f  dx'  dy '  1 1  [x,y  )  j 2- 1  h  (x  -x  ,y  -y  ’)  | 2  •  f  f  dx'  dy  |  h  (x  -x.y  -3 y  )  | 2 


(U) 


Due  to  Schwart’s  inequality.  0£  M{x.y )  <  1  .  Moreover,  using  the  condition  for 
equality  in  Schwartz's  inequality,  the  multiplicative  model  holds,  i.e.  M~ 0  ,  if 
and  only  if  f(x'.y’)  is  constant  over  the  resolution  cell  of  the  system  (where 
/t(x-x'.y-V)  *  0). 

Figures  3  and  4  describe  the  spatial  distributions  of  M(x,y )  and  1 ^  (x.y) 
for  a  sharp  edge  as  well  as  for  slanted  (see  insert  in  Fig .4)  objects,  as  obtained 
from  an  imaging  system  with  unit  magnification  and  a  square  aperture.  It  is 
readily  seen  that  far  from  the  edge  M  approaches  zero.  Also,  as  the  slope  of  the 
object  decreases,  so  does  M  . 

We  have  thus  shown  that  the  multiplicative  model  fails  when  the  object  con¬ 
tains  fine  details  which  cannot  be  resolved  by  the  imaging  system. 


ACKN  OWLEDGEMENTS 


This  work  was  supported  in  part  by  the  Air  Force  Office  of  Scientific 
Research.  M.  Tur  would  like  to  acknowledge  the  support  of  a  Rotschild  and  a  Full- 
bright  Fellowships. 


REFERENCES 


1.  J.  C.  Dainty.  Ed..  Laser  Speckle  and  Related  Phenomena,  Vol.  9  of  Topics  in 
Applied  Physics,  Springer  Verlag  (1975) 

2.  J.  S.  Tjm  and  H.  Nawab,  Opt.  Eng. ,20,  472-480(1981) 

3.  J.  S.  Lee.  Computer  Graphics  and  Image  Processing,  17,  24-32(1981) 

4  A.  K.  Jain  and  C.  R.  Christensen,  Proc.  SPIE  Conf.  on  Applications  of  Speckle 
Phenomena,  Vol.  243,  46-50 (i960) 

5.  H.  Kato  and  J.  W.  Goodman,  Appl.  Opt.,  14,  1813-1823(1975) 

6.  H.  H.  Arsenault  and  G.  April.  J.  Opt.  Soc.  Am..  66,  1160-1163(1976) 

7.  C.  R.  Christensen  ef.  al,  "Object  detectability  in  speckle  noise",  Interna¬ 
tional  Conf.  on  Lasers,  Orlando,  Florida.  Dec.  1978. 

8.  S.  Lowenthal  and  H.  Arsenault,  J.  Opt.  Soc.  Am.,  60,1478-1483(1970) 

J.  W.  Goodman,  in  Ref.  [l]. 


9. 


-7- 


Flgure  Captions 


Figure  1:  The  optical  system.  d(x\y)  and  t(x ,  y)  are,  respectively,  the 
complex  transmittances  of  the  diffuser  and  the  object  transparency. 

Figure  2:  Same  as  Fig.  1  but  with  a  sharp  edge  object.  is  the  part  of 

the  system’s  resolution  cell  (  Rcau  )  which  is  blocked  by  the  object. 

Figure  3:  The  spatial  variation  of  M(x,y)  for  a  sharp  edge  object.  The 
incoherent  image  of  the  sharp  edge  is  also  included.  The  abscissa  is  scaled 
by  the  transverse  dimension  of  Reau  .  X  is  the  wavelength.  /  is  the  focal 
length  of  the  imaging  system  and  o  is  the  size  of  its  square  aperture. 

Figure  4:  Same  as  Fig.  3,  but  for  two  slanted  objects  (see  insert).  While  the 
change  in  object  1  occurs  on  a  scale  comparable  to  Rc,u  .  object  2  is 
fully  resolved  by  the  system.  The  incoherent  impulse  response  of  the  sys¬ 
tem  is  also  included  in  the  insert. 
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